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ABSTRACT 


A basic function in the proper management of repair part inventories 1s the forecasting of 
future demand. The Navy maintains a database of univariate demand data for its repair part 
inventories using a quarterly time interval. Historically, Navy repair part demand forecasting has 
been done using the exponential smpatine procedure. This method 1s a simple and robust means 
of forecasting, however it does not make use of any characteristics of the entire time series such as 
trend, cycles, presence of outliers, or demand clustering. 

This research begins by developing several simple, robust, and dimensionless time series 
features. These features are used to predict the suitability of Box-Jenkins (ARIMA) modeling. 

The ARIMA process is a powerful time series modeling and forecasting technique which possesses 
flexibility for the inclusion of many time series characteristics. This research project develops a 
predictive model of ARIMA suitability using both classical regression and a modern expert-system 
statistical package, ModelQuest. A computationally simple means is presented for determining 
which time series may benefit from the Box-Jenkins methodology. Using ARIMA modeling for 
time series that show significant benefit will provide a more accurate demand forecast and benefit 


inventory management. 
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EXECUTIVE SUMMARY 


This thesis explores a means of determining the suitability of Box-Jenkins, or 
Auto-Regressive Integrated Moving Average (ARIMA), modeling for the demand 
forecasting of Navy repair part data. The ARIMA process is a powerful method of data 
fitting and forecasting that is sensitive to patterns or processes within univariate time 
series. This process is a desirable alternative to other less powerful but computationally 
simpler methods such as exponential smoothing, which is currently in use by the Navy. 
The primary drawback to the use of ARIMA modeling is the increased complexity of the 
model, the intensity of the computation, and the fact that many parts do not benefit from 
its use. It is estimated that 32% of Navy repair parts with recurring demand will show 
significant benefit from ARIMA modeling. 

A method of evaluating which time series will benefit most by ARIMA modeling 
is developed. Several computationally simple features are presented and applied to the 
time series. These features seek to identify the presence of different and relevant 
characteristics within the individual time series. Examples of these features include: 
trend, skewness, seasonality, and autocorrelation. A method of classifying the specific 
time series is generated that takes these features compositely. The suitability of ARIMA 
to a particular demand stream is assessed by processing the set of features through a 
classical regression and an expert-system statistical network program. These methods 
create a fast and simple means of discriminating which time series will benefit most, and 


are therefore worthy of, ARIMA modeling. 
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I. INTRODUCTION 


The Navy Supply System’s primary responsibility is stocking, maintaining, and 
providing the spare parts essential to support naval ships and shore establishments. Proper 
management of these spares 1s necessary to maintain operational effectiveness. Making 
the right choices with regard to which repair parts to stock and at what levels to stock 
them gains importance as the number of ships and their support base decreases. 

A. REPAIR PART FORECASTING 

The Naval Inventory Control Point at Mechanicsburg, PA (NAVICP-Mech) is 
responsible for managing the repair part inventories required by the Navy. Demand 
forecasting is the basic building block of this inventory management. Using historical 
demand data, forecasts are developed for each repair part. These forecasts are then used 
to derive inventory stocking levels. Historically, the method used to generated these 
forecasts has been exponential smoothing. This technique uses a simple weighted average 
between the forecast for the current period and the actual demand for the current period 
to come up with a forecast for the next period. The primary advantages of the exponential 
smoothing process are its simplicity and the low computer data storage requirement. 

Over the past 10 years, however, the aggregate quarterly demand data has been 
kept in a separate database. For each repair part, a single number representing the total 


quarterly demand from all activities was recorded and stored. The aggregate of quarterly 





demand observations for each repair part forms a univariate time series. Typical examples 


of these time series are shown in Figure 1. 
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Figure 1: Examples of Repair Part Demand Histories 


The primary purpose of this database is to capture a demand history, or demand 
stream, for possible future analysis. No covariates of demand have been developed nor 
recorded, however these are certainly worthy of consideration for selected instances in 


future analyses. Possible covariates include: number of ships or number of shore 











part is employed, the criticality of the repair part, the size or handling difficulty of the part, 
the shelf life, and the mean age of the supported equipment. 

The goal of this thesis is to examine and explore the characteristics, or features, of 
these univariate demand streams and develop measures of suitability for Box-Jenkins 
modeling. A means of categorizing the demand process emerges upon studying the 
features of a demand stream. This thesis explores various features with respect to their 
ability to both quantitatively and qualitatively anticipate the suitability of Box-Jenkins, or 
auto-regressive integrated moving average (ARIMA), modeling. 

B. THESIS ORGANIZATION 

Chapter IT provides background regarding exponential smoothing, the time series 
features, and the Box-Jenkins methodology. Chapter III contains the motivating 
methodology of the thesis. It describes the foundations for the research, and presents 
details regarding the calculations performed prior to the analysis. Chapter IV contains the 
analyses, both quantitative and qualitative, and presents classification rules for determining 
ARIMA modeling suitability. Finally, Chapter V provides recommendations and 


conclusions. 





Il. BACKGROUND 


NAVICP-Mech is currently responsible for a line item inventory of over 230,000 
stock numbered repair parts. A large number of these repair parts are inactive and have 
experienced no recent demand, though. For this thesis, a database of the demand histories 
for 139,491 repair parts, comprising those with recent active demand, was extracted for 
study. 

A. EXPONENTIAL SMOOTHING 

Exponential smoothing is a technique that incorporates knowledge of the entire 
history of a process to form a forecast of the next observation. It is commonly used as a 
simple and rational means of time series forecasting. The equation for exponential 
smoothing is 

S,=ay, +(- @)s,, 0<a<l (2.1) 
where S, is the forecast made in period ¢, a is the exponential smoothing coefficient, and 


y, is the demand in period ¢. Recursively substituting into equation (1) results in 
S=ay,+Q-a)ay,,+0-a)yay,,+...+0-a) a y, (2.2) 
This result shows that the forecast for any time t is a linear combination of all previous 


demands observed, with increasing weight given to more recent observations. 











The initial impetus for using exponential smoothing to forecast repair part demand 
was the low requirement of computer data storage and simplicity of calculations. Strict 
boundaries on amount of data archivable in the early computers ruled out other choices. 
B. DEMAND STREAM CHARACTERISTICS 

The shortfall of exponential smoothing is that it does not consider characteristics 
of the demand stream as a whole. There are possibly both long and short term processes 
at work in a time series. For example, features such as trend, cycle, or discrete shifts in 
the process’ mean or variation may be characteristic of a particular demand stream. The 
characteristics of each demand stream can uniquely describe and categorize it. Bissinger 
and Boyarski (1992) present a detailed discussion of the possible patterns in that may 
occur. This thesis develops a basic ability to evaluate the presence of certain distinct 
features. The suitability of alternative modeling methods such as the Box-Jenkins 
methodology can be determined using this ability. 

Cc: BOX-JENKINS FORECASTING 

Two powerful time-series modeling processes combined by Box and Jenkins are 
moving average (MA) and autoregression (AR) (Box and Jenkins, 1970). These 
processes attempt to separate the time series process into its components of an underlying 
process and a white noise process. The white noise process contributes the unpredictable 
random element to the time-series process, while the underlying process can be examined 
and used for analysis and prediction. | 


A moving average process of order q is defined by 





A SAAD Eg it De ee, (2.3) 
where E; is an unobservable error term and the {B,} represent the moving average 
parameters. The interpretation of the MA process of order q is that any given value in the 
series, X;, is directly proportional to the g previous random errors, {E, } , plus some 


current random error E;. This means that the prediction for period ¢ is based only on the 
random errors that have occurred in the g previous periods. 
An autoregressive process of order p is one which recursively satisfies 


Ap eA Xi PAA bt Ag te (2.4) 
where £; is again an error term representing a process with mean of zero and a finite 


variance and the {A,} are the auto-regressive parameters. An autoregressive model of 


order p expresses a time series value as the arithmetic combination of p past series values 
plus a random error term. With p=1, the autoregressive process is identical to an 
exponential smoothing process (Box and Jenkins, 1970). 

Combining these two processes into one with both MA and AR parameters results 
in a model of a stationary autoregressive moving average process (ARMA). A stationary 
process has a constant underlying mean. The ARMA model represents any series value as 
the combination of both past series values and past random error values. The assumption 
of stationarity is frequently an inappropriate one, though, and the regular ARMA process 
may therefore be inadequate. Stationarity can often be induced by using a differenced 


demand history. A single differenced series is attained through the equation: 











Z,=%X,-X,.; (2.5) 
This differencing operation is integrated into the ARMA process by substituting the 
differenced values, Z,, for the original values. The difference process may be repeated as 
often as necessary to provide stationarity. The order of the differencing operator is 
denoted by d. 

The inclusion of the differencing operator to the ARMA process results in an 
autoregressive integrated moving average (ARIMA) model. The ARIMA(7,dq) 
model combines several powerful approaches to time-series analysis into a single package. 
The process is defined by all three parameters: the order of the autoregression, p, the 
order of the integration operator, d, and the order of the moving average process, gq. 

The difficulty in applying the ARIMA process to the NAVICP-Mech demand 
streams is two-fold. First, the time series available for analysis are short ones. The 
ARIMA process works best over long times series that exhibit stationarity, but in the 
repair part database there are a maximum of 40 time intervals available for analysis. 
However, maintaining a more lengthy history will not likely assist the modeling process in 
this case. The database already contains 10 years worth of data — a long period of time 
when considering the total life cycle for a supported system. The addition of demand 
measurements that are over 10 years old will likely provide little insight into the present - 
demand process. 

The second difficulty is the necessary computer time. Since ARIMA will not 


provide significant benefits in every instance, it is prohibitive to calculate directly the 





suitability for all repair parts. Using the statistical software package S-Plus” on a typical 
personal computer system running at 66 Mhz, the calculation of a proper ARIMA model 
to fit a single demand stream requires slightly over 10 seconds of processing time. Using a 
large mainframe computer will certainly reduce the processing time greatly, however the 
time required to process a multitude of demand streams will still be quite large. A better 


way to identify stock numbers which benefit by ARIMA modeling must be found in order 


to take advantage of the modeling benefits it affords. 





Hil. METHODOLOGY 


Since it is infeasible to conduct an analysis of the entire data set provided by 
NAVICP-Mech, we first develop a logical methodology of sampling and evaluation. The 
initial task in defining the methodology is the specification of a sampling plan. The one 
used in this thesis begins by partitioning the data into subsets in which some degree of 
homogeneity can be expected. Two measures are selected for partitioning: the coefficient 
of variation and the number of zeros within the demand stream. Once these cross 
classified subdivisions are established, each one is sampled at random to provide demand 
streams for detailed analysis. 

A measure of effectiveness (MOE) to describe ARIMA modeling suitability is then 
defined and applied to all the sampled demand streams. A means of predicting that same 
effectiveness is developed next. This prediction is based upon the characteristic features 
of each of the demand streams. Each feature is defined based on its ability to measure and 
describe a basic characteristic. 

A. SAMPLING SCHEME 

Not surprisingly, a large portion of the demand streams and the database as a 
whole consists of zeros. This is due to the high percentage of repair part items which are 
maintained as insurance items but which rarely experience demand. Demand streams with 
a high number of zeros are exceedingly difficult to analyze. For the purpose of this thesis, 


a high total number of zeros is viewed as a process with an intermittent random 





component. Only those demand streams with at least 20 positive demands were 
considered for analysis. After the demand streams with high zero counts are filtered out, 
the set of interest, No, comprises 12,000 demand histories. 

The 12,000 observations in N, are viewed in several different ways in search of 
clear partitions. These partitions are chosen to define a sampling scheme, and possibly to 
provide an initial insight into the distribution of stock numbers by feature. The work of 
Boyarski (1995) suggests several different measures to be used in this initial analysis. 
Specifically, the following characteristics are identified as very meaningful: demand mean, 
demand stream length, coefficient of variation, and number of zeros within the demand 
stream. Additional measures that were considered were mean length of non-zero runs, 
mean length of zero runs, median, and the ratio between the mean and the median. 

The distribution in the data for each of these measures was examined and graphed, 
both singly and in bivariate groups, in search of defensible partition values. While most 
characteristics had distributions that displayed a sharp right skewness and no clear 
divisions, the coefficient of variation distribution was unimodal, surprisingly symmetric, 
and centered rather close to one. 

Based on the encouraging distribution of the coefficient of variation, it was chosen 
as the primary criteria for developing a means to partition No. The quantile breakpoints 


for the 12,000 element subset were evaluated and are shown in Table 1. 
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Coefficient of 
Quantile variation value 








Table 1: Empirical quantiles of dataset No 


The five groups resulting from separation by coefficient of variation (c.v.) quantile are 
then examined further. A distinctive feature which provides for a large degree of variation 
within the quantile subsets is the number of zeros within the demand stream, and this 
feature is chosen as a secondary discriminator. 

In using this discriminator, it is again desirable to have roughly equal numbers of 
stock numbers in each bin. Since the spread of the number of zeros varies greatly between 
the different c.v. quantile subsets, secondary breakpoints are chosen in each of the subsets. 
Table 2 shows the final partitioning, divided by range of number of zeros within the time 


series. 


1] 





Subclass: 





Table 2: Divisions of N, by Number of Zeros within the Demand Stream 
An identification scheme for each of the divisions is developed. Each cell is 
identified by two characters. The class, an alphabetic character A through E, corresponds 
to the divisions by coefficient of variation. The subclass, a numeric character | through 5, 
is assigned based on grouping by number of zeros within the class. 
A sample of 45 demand streams was taken from each subset of No. This 
represents roughly a 10% sample from each division, with a total of 1080 demand streams 


sampled. The analysis portion of this thesis relies on these subsets to create and test 


statistical models. 
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B. MEASURE OF EFFECTIVENESS (MOE) 

The base MOE adopted for a time-series model 1s the ratio of the standard 
deviation of the demand stream to the mean squared residuals for the predictive model. 
Algebraically, 


vor (3.1) 


A 
LO. - ys 
n-| 
An effective modeling process is one for which this ratio is above one. 
C; FEATURES 
Several demand stream features are developed and implemented in the course of 
this thesis. Of primary importance while developing these features are the considerations 
of computational complexity, robustness, and dimensionality. A very low computational 
complexity is essential in making the calculation of the features an inexpensive process. 
Robustness implies that the features can both be calculated almost universally and will 
provide meaningful values over all of the different demand streams. Dimensionless 
measures are desirable since their resultant values can be easily compared from demand 
stream to demand stream. When taken compositely, these measures provide for the 


classification of a particular repair part. 
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1. Coefficient of Variation 


The coefficient of variation is defined as x . By measuring the standard 


deviation as a fraction of the mean, the coefficient of variation evaluates how “noisy” a 
process is. Larger values of the ratio correspond to increasingly more random, or noisy, 
processes. 

ze Number of Zeros 

The number of zeros is simply the count of zeros within the time series. In the 
NAVICP-Mech database, there are a large number of demand streams that begin with one 
(or more) zeros. These leading zeros were not considered when counting the number of 
zeros. For the purpose of this thesis, the count of number of zeros begins at the first non- 
zero entry in the time series. 

ae Trend 

Trend is the component of a demand stream that measures the rate of growth or 
decline over time. The feature of trend seeks to quantify the existence and degree of the 
trend component within a series. The development of this feature begins with considering 
the demand series in the form (tf; ,y;) where 7 = 1,...n , t; denotes the specific quarter and ); 
is the demand for the quarter. The demand series is then divided into thirds according to 
time. This means that n = 3k+r where k is a positive integer and r is the remainder term, 


r =0,1,2. Table 3 displays the different divisions possible. 
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Table 3: Series divisions based on remainder 


For the pairs in the lower third, (#; ,y;), where i = 1,...k, let my and y,; denote the 
medians. Similarly, for the pairs in the upper third, (¢; ,.y;), where 7 = n-k+J,...n, let m,3 and 
13 be the respective medians as well. 


From these results, a robust line (McNeil, 1977) is calculated by 


¥3— Yn 


(3.2) 
M3, — My 


In order to convert this to a dimensionless index, the demands are first ordered such that 
Vix = Voy Ses S Yeny (3.3) 
Then, the /00/6 and 500/6 quantiles are estimated. These values are denoted by yz and 


56) Tespectively, and play the role of m,; and m,3 respectively. The dimensionless trend 


index then is given by 


Trend Index = —28— 28 — (3.4) 


¥si6) — Yas6) 
For the trend index, values close to zero support no trend, while values close to +/- 1 


support trend of the appropriate sign. 


15 








4. Peaks 

Peaks are outliers of the demand process. A peak represents a significant and 
abrupt deviation from the established process. A robust measure of peaks was developed 
for this thesis. It utilizes a 10% trimmed mean and variance. Using the trimmed mean and 
variance provides for a magnitude resistant method of measuring outliers. These measures 
are given by 


1 70.9) 





Mon He dM (3.5) 
0.1) 
: ] 79.9) _ : 
£ (0.1) — by YG) —~Y1 (0.1) ) (3.6) 
Nog = I 0.1) 


where n,, represents the number of observations used, and 7,9 1), 0.9) represent the first 


and last subscripts used within the trimmed series 
The measure of peaks is measured in terms of the number of standard deviations 


from the mean. The measure of peaks is made relative to the standardized value 





_ y, Vir 6 
ee ae (3.7) 


i 
5; (0.1) 


By Normal theory, it is expected that about 5% of these measures within a series will be 
larger than 2 and these members are defined as peaks . The feature of interest is the total 


count of the number of peaks within a series. 
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5, Seasonality 


Seasonal variations consist of regular cycles which occur when demand fluctuates 
in a repetitive pattern. The most evident cycle in every day usage is an annual one, and the 
feature of seasonality in this thesis gauges the degree to which repair part demand 
fluctuates according to the four seasons within one year. Other periods of seasonality may 
prove useful, and can be calculated easily through simply modifying the method described 
below. 

A winsorized time series is used in order for this feature to be resistant to outliers. 
The winsorization process replaces the tail order statistics so that outliers do not provide 
excessive influence (Miller, 1986). At the same time, it prone the data loss associated 


with trimming. The winsorized time series is given by 


l _ 1.5%.) —] 


Nosy)? 
Ay = 4 Xs P= No.1) 2--->No9) (3.8) 


nosytt? 2 = Nosy + 1,...,72 


Our feature of seasonality requires that the series be distributed among four 


bins. Let n = 4k +r, r = 0,1,2,3. Then the winsorized series is separated into four 
bins, y,; ,¥2; »¥3,; and y,;, with each containing & or k+1 entries. A series element is a 
member of bin y, ; when it can be written as y;,;+4/;-1) for 7 = 1,2,3,4 andj = 1,2,...k. 

For each of the bins, ¥, and 5, is calculated. In effect, a one-way analysis of variance 


(ANOVA) is performed to produce the within and total sums of squares. The coefficient 


of determination is given by 
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Re (3.9) 
SumOfSQuar ess orc 


This serves as the measure of seasonality, where values close to one indicate high 
seasonality, while values close to zero indicate no seasonality. 

6. Runs 

The feature of runs measures the tendency within a demand stream to have 
consecutive observations either above or below the median of the overall process. A run 
is defined as a collection of successive observations all on the same side of the process 
median. This measure of runs is useful in representing the amount of clustering of 
observations that occurs. 

The clustering property is evidenced by groupings, or clumps, of observations 
interspersed by periods of no or little demand (Bissinger and Boyarski, 1992). Small 
values of the runs feature correspond to a tendency to cluster while large values indicate a 
more oscillatory process. This feature should always be considered in conjunction with 
the trend feature, since a large trend component may lead to a false indication of no 
clustering. 

In measuring runs, first calculate the median of the demand stream. If any demand 
observations are equal to the median, they are discarded for the computation. The number 
of runs above and below the median are then counted. These are denoted by R; and Ra, 


respectively. The sum of these two measures is denoted by R. Two additional necessary 
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measures are 72, the total number of observations above the median, and 72, the total 
number of observations below the median. 
According to Gibbons (1992), the expected value and variance of R are given as 


E[R] = 1+ ok (3.10) 


> 
n, + Hn, 


2n, n,(2n, n, ~ nm ~N)) 


Var(kR) = 3.11 
(1, +n,) (1, +N, — 1) 
The feature of runs can then be standardized into a dimensionless measure 
ee sae (3.12) 


J Var(R) 


Te Skewness 
Skewness is the tendency of a distribution to be asymmetric. To evaluate 


skewness, the comparison is made between the mean and the median of a series. Simply, 


I, = i (3.13) 


where y is the arithmetic average demand and m, is the median demand. For J, values 


an index of skewness is given by 


greater than one, the process is skewed to right. Conversely, values less than one indicate 


a process skewed to left. 
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8. Autocorrelation 

Autocorrelation is a measure of how strongly time series values are correlated to 
each other over time. A nonparametric measure of the overall autocorrelation based on a 
single lag element is chosen. The measure presented by Gibbons (1992) is used. 

First, assign each observation y; a rank 7; relative to all the other elements of the 
series. That is to say, the lowest ; is assigned a rank of one and the largest is assigned a 


rank of n. Then form the following measures: 


— - 
RVN = =—Gite (3.15) 
Dale - = 


i=] 
If there are no ties within the ranks, this simplifies as 


12 NM 


RVN = mn? — 1) (3.16) 
The expected value of this measure is 2, and its variance is given by 
am 4(n — 2)(Sn* - 2n - 9) (3.17) 
5n(n + 1)(n — 1)° 
This allows RVN to be normalized according to 
RVN om = RIN = (3.18) 
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This result is the feature of autocorrelation, with small values indicative of positive 


autocorrelation and large values indicating negative autocorrelation. 
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IV. ANALYSIS 


The characteristic features are first evaluated and stored for each of the 1080 
elements of the dataset. The program to calculate the features and summarize them into a 
convenient database was written in Turbo-Pascal. When compiled, the Turbo-Pascal code 


is extremely efficient in terms of speed of calculation. Less than 30 seconds of total 


the entire dataset. 

The MOE of ARIMA suitability is calculated using the S-Plus” statistical package. 
The MOEs for exponential smoothing as well as for three different orders of ARIMA 
models were calculated. Due to the large number of different ARIMA models available, 
the three models were chosen to provide as great a range as possible. The three models 
chosen are: a first order (1,1,1) model, a second order (2,2,2) model, and a third order 
(3,2,3) model. In the third order model, the order of integration 1s left as two since 
integrations of order three (or higher) are almost never necessary (Hoff, 1983). The 


results are shown in Figure 2. 


computational time are required to calculate and store the values for all of the features for 
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Model type 


Figure 2: Boxplots of MOE Values for Four Different Modeling Techniques 


Figure 2 shows that there are a number of demand histories for which ARIMA 
modeling provides an appreciable benefit over exponential smoothing. At the same time, 
however, there are a number of demand histories for which exponential smoothing 
provides nearly as much predictive power as any of the other models. The mean MOE 
improvement when compared to exponential smoothing is: 0.15 for the ARIMA(I,1,1) 


model. 0.28 for the ARIMA(2,2,2) model, and 0.41 for the ARIMA(3,2,3) model. 
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While it would be desirable to have MOE improvement in separate types of time 
series using separate levels of ARIMA modeling, such does not appear to be possible. 
Comparison of the MOE results shows that there was one general group of time series 
which benefited from ARIMA modeling. Where ARIMA(1,1,1) provides an appreciable 
benefit over exponential smoothing, higher order ARIMA models generally provide an 
additional benefit over the single order model. In those cases where ARIMA(1,1,1) 
provides little benefit, the other models usually provide little improvement as well. This 
observation is evidenced again in the statistical network analysis. 

A. STAR PLOTS 

A better understanding as to which features of a demand history are indicative of 
ARIMA model suitability is desired. An initial means of exploring this is through star 
plots as described by Chambers, Cleveland, Kleiner, and Tukey (1983). A star plot 
converts multi-attribute data into graphical stars. The lengths of each specific ray of the 
star corresponds to the magnitude of that attribute relative to the same attribute within the 
other stars. A pattern may emerge of features that are indicative of model suitability by 
grouping star plots by MOE level and comparing those demand histories that benefited 
from ARIMA modeling. A diagram of the star plot format and examples of star plots with 
likely ARIMA modeling benefit are shown in Appendix A. 

Examination of the star plots shows that several of the features appear to have 
ARIMA predictive power. A large measure of autocorrelation paired with an indication of 


trend is a common element in many of the demand histories that benefited most from 
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ARIMA. Coefficient of variation does not seem to play a significant role. A moderate 
measure of skewness was usually evident in the histories with larger MOEs. Several of 
the star plots were difficult to visually identify as ARIMA candidates, due to the lack ofa 
definitive common element. There is also a wide disparity between the star plots of those 
histories which experience a moderate improvement in the MOE and those with little 
improvement. 

The inclusion of additional features in the star plots was explored on the computer. 
The features of variance, mean, mean length of zero runs, and total number of 
observations within the demand stream were added. None of these additional features 
provide any improvement to the visual categorization of the star plots, and often detracted 
from the categorization. The additional features were subsequently removed from further 
consideration. 

An interesting result from the examination of the star plots is that the samples 
taken from the E class, representing those with a coefficient of variation greater than 1.5, 
have little similarity with those from the smaller coefficient of variation subsets. Elements 
from the E class with a high MOE measure have a variety of star plot shapes, and are 
generally very different from those with a high MOE from the other subsets. This suggests 
that there may be an undiscovered process at work in the high coefficient of variation 


demand streams that is obscuring ARIMA model suitability. 
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B. CLASSICAL REGRESSION 
The standard method for determining the relationship between predictor variables, 
namely the demand features, and resultant variable, the MOE, is through classical 
regression. Regression analysis is conducted using the eight explanatory variables from 
the star plots: coefficient of variation, number of zeros, trend, peaks, seasonality, runs, 
skewness, and autocorrelation. The first regression tested is a simple additive model of 
the explanatory variables, as shown in the Equation 4.1. Initially, the MOE values for 
_ ARIMA(1,1,1) models are used as a resultant measure. 
y = By + Bix, + BX, + Bax, + Bax, + BsxX5 + BeX— + B>X; (4.1) 
where y is the MOE measure, 

x, 18 the coefficient of variation, 

X2 is the number of zeros within the series, 

x3 is the trend feature, 

x4 is the measure of peaks, 

x5 1s the seasonality feature, 

Xs 18 the measure of runs, 

X7 1s the skewness feature, 

Xg is the autocorrelation measure, 

and /; are the regression coefficients 

The regression results for the data set predicting ARIMA(1,1,1) suitability are very 

modest. The resultant p-value of 0.0 indicates strong evidence of a relationship between 
the predictors and the resultant, but the R’ value of 0.16 implies that only a small fraction 
of the MOE can be successfully predicted. 


The post-regression analysis provides a likely reason for this poor performance. 


Examination of the raw residuals shows several that are outlier candidates. After 
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calculating their Studentized residuals and Cook’s distances, it is clear that these are 
highly influential outlier points. Interestingly, the outliers are primarily grouped together 
within the E class. Further investigation indicates that almost all the elements of the E 
class had significantly larger Cook’s distances compared to the A through D classes. 
Based on this result, and the results from the qualitative analysis of the star plots which 
show large variability within the E class, a structurally identical regression is performed 
using classes A through D only. 

The results of the second additive regression are very promising. The p-value 
remains at 0.0, while the R’ value increases to 0.51. The post-regression analysis reveals a 
lingering bow in the residuals vs. fitted values plot, but otherwise no surprising results. 
Excluding the class with the largest coefficient of variation results in a significant 
improvement in predictive power. It is evident that the E class will require special 
treatment, since it behaves noticeably differently than the other subsets. For the purposes 
of this thesis, modeling will be done using only the classes A through D. 

Consideration was given as to whether the response variable, MOE measure, 
should be transformed to provide a better fit. The primary indication that a better fit might 
be obtained is the bow in the residuals plot. We make use of the quantitative method Box 
and Cox (1964) developed to find the optimal transformation. After fitting a full model to 
the data, the Box-Cox transformation algorithm is applied to the model to obtain a 95% 
confidence interval for the optimal transformation. Based on this result, a power 


transformation with A= -3.5 is applied to the full model. This results in a p-value of 0.0 
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and an R? value of 0.65. The R’ values between the transformed and untransformed 
models cannot be directly compared, since they model different responses. However, the 
results of the transformed regression can be “untransformed” and a comparative R’ 
generated using the ratio of the within sum of squares to total sum of squares. This 
calculation results in an R’ of 0.68, representing an average variance reduction from the 
original R’ of 0.51. The post-regression analysis shows a much improved distribution of 
residuals with the transformed model as well. 

A determination as to which, if any, of the predictor variables are extraneous in the 
model is appropriate in validating the model. A stepwise regression process is performed. 
This thesis uses Efroymson’s stepwise method, which is similar to forward selection. The 
significant difference with forward selection is that after each predictor addition, partial 
correlations are considered to see if any of the variables in the subset should now be 
dropped. By either adding or subtracting terms at each step, an efficient model with a 
minimal number of terms results. All the single variable terms as well as the second order 
interaction terms between predictors are included in the stepwise regression. The 
interaction terms are significant in that they consider the relationships between predictive 
variables. From the stepwise regression, the final regression model for predicting single 


order ARIMA coefficients 1s 


-3.5 
y i Bot BX Poke +P ke tpt + Ps + Pass + 


(4.2) 
Bj X_X, + ByX,X_ + PoX4X, + ByXsXq 


where y is the MOE measure, 
x, is the coefficient of variation, 
X2 1s the number of zeros within the series, 
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x3 is the trend feature, 

x4 is the measure of peaks, 

x5 is the seasonality feature, 

X¢ is the measure of runs, 

x7 1s the skewness feature, 

Xe is the autocorrelation measure, 

and / are the regression coefficients. 
Noteworthy in these results is the fact that all of the features previously developed 
contributed to the overall regression equation, either singly or as a part of an interaction 
term. The untransformed R’ value for the stepwise regression result is 0.7734. The 
detailed results of the regression for the ARIMA(1,1,1) resultant are shown in Appendix 
B. 

Graphing the tranformation inverted fitted values vs. their residuals shows an 
interesting result. Resultant values at the high-end of the MOE scale, above approximately 
1.6, are predicted to be significantly higher than they actually are. This result is due to the 
large power of the transformation, which skews the values when the transform in inverted. 
With this knowledge, actual predictions from the regression model would have an 
effective R’ greater than the evaluated untransformed one of 0.7734 since a compensation 
can be made to larger values of predicted MOE. The plots of transformation inverted 
fitted values vs. residuals are shown in Appendix C for all three cases that were ARIMA 
modeled. 

A similar regression analysis as with the single order ARIMA is conducted using 


the ARIMA(2,2,2) and ARIMA(3,2,3) MOE values and calculated features. The method 


of analysis is the same as for the ARIMA(1,1,1) case. In brief, the simple additive model 
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with an untransformed result is first developed. A transformation for the result variable, 
using the Box-Cox algorithm, is then considered. Finally, a stepwise regression is 
conducted using the transformed result variable and all the predictor variables as well as 
their interaction terms as explanatory variables. At all steps in the regression process, the 
results are critically evaluated for influential points, outliers, or other regression anomalies. 


For the ARIMA(2,2,2) case, the final equation is of the form 


yr = Po + Pixs F Paks + Pyx¢. FP Xe FP Xpks +P XXet 
B7X4X_ + BpX5X5 + PoxX5X, 


(4.3) 
The regression equation for the ARIMA(2,2,2) case results in an R’ of 0.6908 when the 
transformation is inverted. Inspection of the plot of untransformed fitted values vs. 
untransformed residuals shows a similar result as with the ARIMA(1,1,1) case. There are 
three points at the high-end of the MOE values which contribute the bulk of the error in 
prediction. As a result of the powerful transformation in Equation 4.3, these values are 
predicted to be significantly higher than their actual values. The plot of untransformed 
fitted values vs. untransformed residuals is shown in Appendix C. The detailed results of 
the regression model for the ARIMA(2,2,2) case are shown in Appendix B. 


For the ARIMA(3,2,3) case, the results are again similar to the previous cases. 


The resulting regression equation is of the form 


yo" = B, + B,x,x, + B,x,x, + B,X,x, + Bux, xX, + Byx.X_ + BxX,X, + 


(4.4) 
PahaheF PeXak yt Pokake © Pere Pikes Pike 
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previous two models, this R’ is greatly affected by a few demand streams with large MOE 
values which are overestimated by the model. The decreased power of transformation 
used in the ARIMA(3,2,3) model results in a noticeably lessened distortion effect in the 
tails of the predicted values. The specific regression details for the ARIMA(3,2,3) model 
are shown in Appendix B. 

All three of the regression models are highly predictive of their respective MOE 
measures. However, there is a general tendency to overpredict MOE measures when the 
actual MOE measure is above roughly 1.5. This should not be a significant problem in 
predicting ARIMA suitability, though, since any MOE value above 1.5 should be highly 
suggestive of ARIMA benefit to the time series. 

There are, however, other shortcomings to regression analysis. A resulting 
regression model can be easily applied — but the development of the model is a detailed 
and complicated process. Due to the highly specialized nature of the model, an entirely 
new analysis would be required to consider suitability of other ARIMA models. A 
possible alternative to regression analysis is the use of statistical network software. 

C. STATISTICAL NETWORK 

This thesis uses the expert system statistical network software package 
ModelQuest™. ModelQuest™ develops polynomial-order networks of predictive 
variables to model a result variable (ModelQuest, 1995). In many aspects, ModelQuest™ 


performs a similar analysis to regression, however the model building process is entirely 
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performs a similar analysis to regression, however the model building process is entirely 
automated. The goal of using ModelQuest™ as an alternative to regression is to find an 
easier and faster method of generating new models, while maintaining predictive power. 
ModelQuest™ determines both the network structure as well as the coefficients 
within the network. This is in contrast to regression, for which the network structure 
mie be pre-defined. The additional range of freedom that ModelQuest™ has in selecting 
network structure is managed through a predicted squared error (PSE) equation (Barron, 
1984). The PSE for a model is the sum of the fitted squared error (FSE) and a complexity 
penalty (KP). ModelQuest™ defines the optimal solution for a given situation as that one 


which minimizes the PSE value. This is shown schematically in Figure 3. 


Model Performance Error 


Best Network PSE 





Increasing Model Complexity 


Figure 3: Relationship Between PSE, FSE, KP, and the Optimal Solution 


FSE is comparable to the residual sum of squares for a particular regression model. 


ModelQuest™ defines KP in the following manner, 
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KP = (CPMY(=Y(s,") (4.5) 


where CPM is the complexity penalty multiplier 
K is the number of coefficients within the network 
Nis the number of data elements being modeled 
s,’ is an estimate of the model error variance 
The value of the complexity multiplier (CPM) is user-definable. By adjusting the CPM, 
the complexity of a statistical network is controlled — an increased CPM results in a 
decreased complexity model. 

In practice, managing CPM is essential to ensuring that the resulting model is 
neither over- or under-specific. This is done by comparing the performance of the model 
between two different datasets. The training dataset is used to develop the model, while 
the evaluation dataset is used to evaluate the developed model. If the performance of the 
model is significantly different between the training and evaluation datasets, then the 
model has overfit the training data. For this thesis, 75% of the entire database of features 
and MOE values were set aside for the training dataset and the remaining 25% formed the 
evaluation database. The selection of the particular elements for the training and 
evaluation sets was done randomly by ModelQuest™. 

ModelQuest™ creates a layered network model. Polynomial relationships are 
created between variables within each layer of the network. These relationships may 


include the addition of constant terms, as well as polynomials of order up to three. A 


schematic network is shown in Figure 4. 
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Figure 4: Example Three-Layer ModelQuest™ Network 


At each “box” operator of the network in Figure 4, a polynomial of order three (or less) 
combines the left-side input variables of that box to form a result. Where an input variable 
is the output from a previous box, that previous polynomial is treated as if it were a single 
variable. In other words, the order of the polynomial equation which represents the entire 
network may be three times the size of n, where n is the total number of layers in the 
network. ModelQuest™ allows the definition of a strict upper bound for the number of 
layers within the network as well as the number of terms within the first layer of the 
network. 

Other than CPM and the dataset itself, ModelQuest™ requires no other 
information or intervention to develop the network. The resulting networks are more 
complicated than regression results, as will be seen, but they have the benefit of 
completely automated formation. For the statistical network models in this thesis, several 
different values of CPM were tried for each network. The best network model is arrived 
at by comparing the performance within a model for the different datasets as well as 


comparing the performances between models with different CPM. 
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The structure of the statistical network model for the ARIMA(1,1,1) case is shown 


in Figure 5. 


Number of Zeros 
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Autocorrelation 


ARIMA(1,1,1] MOE 
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Runs 


Figure 5: Statistical Network for ARIMA(1,1,1) Case 


In comparison to the regression model, the statistical network uses only five of the 
eight features to predict MOE in this case. This is compensated for by the additional 
complexity of the sixth-order polynomial which results from two triplet nodes. The R’ 
which results from this model is 0.8054 for the training dataset and 0.78032 for the 
evaluation dataset. Details of the statistical network as well as the diagnostics for this 
model are shown in Appendix D. 

The statistical network for the ARIMA(2,2,2) case is more complex than the single 
order ARIMA model. This model uses six of the eight predictive features, adding the 
feature of seasonality to the five features used in the ARIMA(1,1,1) network. The 


structure of the model is shown in Figure 6. 
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Figure 6: Statistical Network for ARIMA(2,2,2) Case 


The structure of the ARIMA(2,2,2) network is very similar to that of the 
ARIMA(1,1,1) case. The primary difference between the structure of the two is the 
addition of a third triplet operation in the ARIMA(2,2,2) case. For the second order 
ARIMA case, the R’ value is 0.5861 for the training set and 0.5869 for the evaluation set. 
The detailed form of the statistical network and the diagnostics of the model are shown in 
Appendix D. 

Interestingly, the structure of the statistical network for the ARIMA(3,2,3) case is 
identical to the ARIMA(1,1,1) case. The only difference between the two networks are 
the coefficients of the predictors. For the third order ARIMA model, an R’ of 0.5751 for 
the training set and 0.6131 for the evaluation set results. The detailed form of the network 
along with the diagnostics of the third order model are shown in Appendix D. 

Applying the network structure of Figure 6 to the second order ARIMA case gives 


an R’ of 0.5097. This provides strong evidence that there is a positive correlation between 
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MOE improvement among the ARIMA processes. Also, the structure shown in Figure 6 
seems to be a good general structure for modeling Box-Jenkins suitability. 

Both the ModelQuest™ and regression models exhibit a high predictive power of 
ARIMA MOE given the developed features. The choice as to which to use in actual 
application depends on the situation. ModelQuest™ requires only a minimum of training 
and a standard personal computer, however the results are not as detailed as those 
available in post-regression analysis. The regression process, in comparison, requires 
careful analysis at each step in the process. The results from regression are extremely 
detailed. However, both regression and ModelQuest™ provide excellent results when 
properly applied. 

D. ARIMA FORECASTING 

It is difficult to develop a meaningful measure of the improvement provided by 
ARIMA modeling over exponential smoothing in actual forecasting. The large variety of 
time series present in the NAVICP database complicate forecast comparisons. Some 
simple forecasts are made to quantify the ARIMA benefit. The first n-/ quarters of data 
are used to forecast the n™ quarter using both exponential smoothing and the ARIMA 
models. This is done for those time series from the analysis dataset which have an 
ARIMA MOE coefficient equal or greater than 1.25. The statistic that is measured to 


quantify improvement is: 
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where F* is the ARIMA forecast for period n, F.” is the exponential smoothing forecast 


for period n, and X, is the actual demand for period n 

This ratio measures the improvement in forecasting from a Box-Jenkins model 
compared to exponential smoothing. The median improvement measure is: 1.10 for the 
ARIMA(1,1,1) model, 1.17 for the ARIMA(2,2,2) model, and 1.22 for the ARIMA(3,2,3) 
model. These results indicate that there is improvement to be gained from forecasting 
with Box-Jenkins models. 

Forecasting from ARIMA models is noteworthy in that predicted demands of zero 
are possible. Exponential smoothing is incapable of producing a forecast of zero, however 
the ARIMA models produce these forecasts recurringly for the NAVICP data. These zero 
forecasts validate to be important due to the large number of repair parts characterized by 
frequent zero demand. 

It is estimated that roughly 8% of NAVICP stock numbered parts would benefit 
from ARIMA modeling. Although this is a seemingly small fraction, 1t represents a more 
sizeable percentage of those repair parts which have recurring sosive demand. These are 
the items which make up the vast majority of all demands, and will correspondingly 


provide a larger overall benefit with improved forecasts. 
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V. CONCLUSION 


This thesis developed eight features to assist in classifying repair part demand 
streams. The features are notable in their ease of computation, robustness, and lack of 
dimension. When taken compositely, the features provide a means of predicting Box- 
Jenkins modeling suitability. 

A. RESULTS 

Both classical regression and expert system software were used to relate the time 
series features to the ARIMA suitability MOE. While arriving at different networks by 
different means, these two methods both result in a high predictability of Box-Jenkins 
suitability. By presenting two different approaches, there is versatility in implementation 
of the ARIMA suitability methodology. The repair parts which will likely benefit from 
Box-Jenkins modeling can be quickly identified through the models developed within the 
thesis. 

B. AREAS FOR FURTHER STUDY 

It would be beneficial to quantify the exact benefits of ARIMA modeling and 
forecasting for Navy repair part data. Such a study would possibly include a look into 
which types of repair parts receive benefit from Box-Jenkins and what the result of an 
improved forecasting technique is. Intuitively, an improvement in forecast accuracy leads 


to decreased risk of a stock out which will allow the safety stock level to be decreased. 


4] 


Further research into demand stream characteristics may also provide interesting 
results. This thesis has developed several features, however additional features likely exist 
as well. A study of the processes within the high coefficient of variation subset would 
directly add to the work of this thesis. For the entire dataset as well, additional features 
would enhance both classification and ARIMA predictive power. 

Sensitivity analysis of the ARIMA models may also be conducted as part of a 
comparison of the efficacy of different models. The versatility of Box-Jenkins through the 
different orders of autoregression, integration, and moving average as applied to the 
NAVICP data is unexplored. There are certainly tradeoffs between different models, and 


research into what they are and how significant a role they play would prove beneficial. 
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APPENDIX A. SAMPLE STAR PLOT DIAGRAMS 
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APPENDIX B. REGRESSION RESULTS AND PLOTS 
ARIMA(1,1,1) case: 


Formula: X9 * (-3.5) ~ X3 + X5 + X6 + X1*X3 + X2*X3 + X2*X5 + X2*X8 + X3*K8 


X4*X7 + X5*X6 
Residuals: 
Min 1Q Median 3Q Max 
-0.5841 -0.09582 0.01352 0.1115 0.4589 
Coefficients: 
Value Std. Error tvalue  Pr(>t)) 
(Intercept) 0.6096 0.0383 15.9132 0.0000 
x3 -0.3321 0.0747 -44476 0.0000 
X5 1.7955 0.4373 4.1063 0.0000 
x6 0.0866 0.0064 13.4368 0.0000 
X1:X3 0.1534 0.0652 2.3525 0.0189 
X2:X3 0.0154 0.0036 4.3194 0.0000 
X2:X5  ~=0.0581 0.0159 3.6501 0.0003 
X2:X8 -0.0012 0.0005 -2.3047 0.0214 
X3:X8 = =-0.1059 0.0124 -8.5294 0.0000 
X4:X7 ~—--0.0043 0.0013 -3.4195 0.0007 
X5:X6 ~~ -0.3533 0.0765 ~4.6191 0.0000 


Residual standard error: 0.1689 on 867 degrees of freedom 
Multiple R-Squared: 0.6848 
F-statistic: 189.5 on 10 and 872 degrees of freedom, the p-value is 0 
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ARIMA(2,2,2) case: 


Formula: X9 “ (-3.0) ~ X3 + X5 + X6 + X8 + X2*X3 + X3*X8 + K4*X7 + K5*X6 + 


3Q 


t value 
11.4773 
~3.0376 
2.3565 
12.1485 
-1.6056 
6.4248 
-4.6574 
-3.7216 
-4.7656 
3.4904 


xX5*X8 
Residuals: 
Min 1Q Median 
-0.5093 -0.09787 0.003696 0.117 
Coefficients: 
Value Std. Error 
(Intercept) 0.5414 0.0472 
V3 -0.1713 0.0564 
V5 1.2312 0.5225 
V6 0.0844 0.0069 
V8 -0.0196 0.0122 
V2...V3 0.0120 0.0019 
V3...V8  -0.0859 0.0185 
V4...V6  -0.0017 0.0004 
V5...V6 = -0.3762 0.0789 
V5...V8 0.3168 0.0908 


Max 
0.4963 


Pr(>|t!) 
0.0000 
0.0025 
0.0187 
0.0000 
0.1087 
0.0000 
0.0000 
0.0002 
0.0000 
0.0005 


Residual standard error: 0.1676 on 864 degrees of freedom 


Multiple R-Squared: 0.6071 


F-statistic: 148.3 on 9 and 864 degrees of freedom, the p-value is 0 
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ARIMA(3,2,3) case: 


Formula: X9 * (-2.4) ~ X1*X2 + X1*X8 + K2*X3 + K2*XK4 + X2*X6 + X3*X4 + 
X3*X6 + X3*X7 + X3*X8 + K5*X6 + K6*X7 + XO*X8 





Residuals: 
Min 1Q Median 3Q Max 
-0.6269 -0.1018 0.0113 0.1024 0.4477 

Coefficients: 
Value Std. Error tvalue  Pr(>t|) 
(Intercept) 0.6457 0.0200 32.2421 0.0000 
V1..V2 0.0052 0.0039 1.3483 0.1779 
V1...V8  -0.0174 0.0145 -1.1954 0.2323 
V2...V3 0.0073 0.0036 2.0326 0.0424 
V2...V4  -0.0004 0.0004 -1.0806 0.2802 
V2...V6  -0.0002 0.0007 -0.2698 0.7873 
V3...V4  -0.0087 0.0046 ~1.8947 0.0585 
V3...V6 0.0159 0.0094 1.6863 0.0921 
V3...V7 = -0.0240 0.0183 -1.3160 0.1885 
V3...V8  -0.1124 0.0133 -8.4247 0.0000 
V5...V6 0.0212 0.0184 1.1502 0.2504 
V6...V7 0.0042 0.0030 1.4131 0.1580 
V6...V8 0.0085 0.0026 3.2751 0.0011 


Residual standard error: 0.1592 on 866 degrees of freedom 
Multiple R-Squared: 0.5337 


F-statistic: 82.59 on 12 and 866 degrees of freedom, the p-value is 0 
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APPENDIX C. UNTRANSFORMED REGRESSION RESIDUAL 
PLOTS 


ARIMA(1,1,1) Regression 


Transform reversed residuals 





Transform reversed fitted 


ARIMA(2,2,2) Regression 


Transform reversed residuals 





Transform reversed fitted 


5] 


Transform reversed residuals 


ARIMA(3,2,3) Regression 





2. 3 


Transform reversed fitted 
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APPENDIX D. STATISTICAL NETWORK RESULTS 


Summary Statistics for output ‘ARIMA(1,1,1) MOE' 


evaluation 
number of observations 22] 
average absolute error 0.070402 
absolute error standard deviation 0.067067 
average squared error 0.0094340 
normalized mean squared error 0.22933 
squared error standard deviation 0.019915 
maximum absolute error 0.41618 
database output minimum 0.91301 
database output maximum 2.2221 
database output mean 1.1318 
database output standard deviation 0.20327 
network output mean 1.1344 
network output standard deviation 0.19955 
R-squared 0.78032 


root of predicted squared error 
predicted squared error 
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trainin 


662 


0.013497 


0.91883 
4.2954 
1.1402 — 
0.26278 


0.15197 
0.023094 





Summary Statistics for output 'ARIMA(2,2,2) MOE' 


evaluation trainin 
number of observations 218 656 
average absolute error 0.13389 
absolute error standard deviation 0.15605 
average squared error 0.042165 0.025480 
normalized mean squared error 0.34533 
squared error standard deviation 0.11735 
maximum absolute error 1.0232 
database output minimum 0.93120 0.91780 
database output maximum 4.1915 4.7695 
database output mean 1.2384 12273 
database output standard deviation 0.34928 0.30813 
network output mean 1.2127 
network output standard deviation 0.26087 
R-squared 0.65635 
root of predicted squared error 0.20524 
predicted squared error 0.042124 
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Summary Statistics for output 'ARIMA(3,2,3) MOE' 


evaluation 
number of observations 220 
average absolute error 0.18444 
absolute error standard deviation 0.24832 
average squared error 0.095399 
normalized mean squared error 0.42826 
squared error standard deviation 0.39671 
maximum absolute error 2.1067 
database output minimum 0.96086 
database output maximum 5.2370 
database output mean 1.3583 
database output standard deviation 0.47302 
network output mean 1.3531 
network output standard deviation 0.27410 
R-squared 0.61307 


root of predicted squared error 
predicted squared error 
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0.054653 


0.97399 
3.958] 
1.3359 
0.35894 


0.26048 
0.067850 
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